A novel method to calculate compliance and airway resistance in ventilated patients

Background The respiratory system’s static compliance (Crs) and airway resistance (Rrs) are measured during an end-inspiratory hold on volume-controlled ventilation (static method). A numerical algorithm is presented to calculate Crs and Rrs during volume-controlled ventilation on a breath-by-breath basis not requiring an end-inspiratory hold (dynamic method). Methods The dynamic method combines a numerical solution of the equation of motion of the respiratory system with frequency analysis of airway signals. The method was validated experimentally with a one-liter test lung using 300 mL and 400 mL tidal volumes. It also was validated clinically using airway signals sampled at 32.25 Hz stored in a historical database as 131.1-s-long epochs. There were 15 patients in the database having epochs on volume-controlled ventilation with breaths displaying end-inspiratory holds. This allowed for the reliable calculation of paired Crs and Rrs values using both static and dynamic methods. Epoch mean values for Crs and Rrs were assessed by both methods and compared in aggregate form and individually for each patient in the study with Pearson’s R2 and Bland–Altman analysis. Figures are shown as median[IQR]. Results Experimental method differences in 880 simulated breaths were 0.3[0.2,0.4] mL·cmH2O−1 for Crs and 0[− 0.2,0.2] cmH2O·s· L−1 for Rrs. Clinical testing included 78,371 breaths found in 3174 epochs meeting criteria with 24[21,30] breaths per epoch. For the aggregate data, Pearson’s R2 were 0.99 and 0.94 for Crs and Rrs, respectively. Bias ± 95% limits of agreement (LOA) were 0.2 ± 1.6 mL·cmH2O−1 for Crs and − 0.2 ± 1.5 cmH2O·s· L−1 for Rrs. Bias ± LOA median values for individual patients were 0.6[− 0.2, 1.4] ± 0.9[0.8, 1.2] mL·cmH2O−1 for Crs and − 0.1[− 0.3, 0.2] ± 0.8[0.5, 1.2] cmH2O·s· L−1 for Rrs. Discussion Experimental and clinical testing produced equivalent paired measurements of Crs and Rrs by the dynamic and static methods under the conditions tested. Conclusions These findings support to the possibility of using the dynamic method in continuously monitoring respiratory system mechanics in patients on ventilatory support with volume-controlled ventilation. Supplementary Information The online version contains supplementary material available at 10.1186/s40635-022-00483-2.

Page 2 of 14 Gutierrez Intensive Care Medicine Experimental (2022) 10:55 The respiratory system (rs) static compliance (C rs ) and airway resistance (R rs ) are calculated during volume-controlled (VC) mechanical ventilation with a breath-hold maneuver at the end of quiet inspiration (static method) [1]. Under these conditions C rs = V tidal /(P plateau -PEEP a ), where V tidal = tidal volume, P plateau = breath-hold P aw , and PEEP a = applied positive end expiratory pressure. Similarly, R rs = (P peak -P plateau )/F aw , where P peak = peak inspiratory pressure and F aw is airway flow measured just prior to breath-holding [2]. A reliable method to calculate C rs and R rs automatically, without the need of an inspiratory hold, would have great utility in monitoring the adequacy of ventilatory support. One approach previously tried is the multiple least squares fit (LSF) technique [3,4], where measures of P aw , F aw , and lung volume change (ΔV) are fitted to the equation of motion of the respiratory system. Another is the expiratory time constant τe method [5] where equations for C rs and R rs are developed assuming mono-exponential lung volume release [6]. Both methods require the use of complex computational techniques and absent respiratory muscle effort.
Described is a method to calculate C rs and R rs during insufflation in the presence of airflow (dynamic method) that combines frequency analysis of the airway signals with a novel numerical solution of the equation of motion. The method was validated experimentally with a one-liter test lung. It was also validated clinically using previously acquired F aw and P aw signal data from patients on VC ventilation displaying end-inspiratory holds. This allowed for the reliable calculation of paired C rs and R rs values using both static and dynamic methods.
Theoretical development. The time-dependent equation of motion of the respiratory system is: This equation, based on the one-compartment model of Otis et al. [7], assumes constant values for C rs and R rs . The measured airway pressure P aw (t) represents the sum of the ventilator and respiratory muscles applied pressures P vent (t) and P mus (t), respectively. Opposing them are the elastic, resistive, and inertial components of the respiratory system. V(t) represents the time-dependent lung volume; ΔV(t) is the insufflation lung volume at time t, equal to t 0 F aw (t)dt ; I is the respiratory system inertia; and PEEP i the intrinsic PEEP [8].
Assuming passive insufflation (P mus = 0), negligible PEEP I, and ignoring the effect of the inertia term [9] , Eq. (1) becomes: It is possible to solve numerically this indeterminate equation with two unknowns, C rs and R rs , by first developing a solution matrix for each set of P aw (t k ), ΔV(t k ), F aw (t k ), and PEEP a values measured at successive times t k during insufflation. The elements of the solution matrix are calculated by substituting the measured values for ΔV(t k ), F aw (t k ), and PEEP a into Eq. 2 and alternately applying a range of physiologically plausible values for C rs (C 1 … C n ) and R rs (R 1 … R n ).
For example, applying a range of (C 1 … C n ) values from 10 to 100 mL·cmH 2 O −1 and 1.0 to 50.0 cmH 2 O·s·L −1 for (R 1 … R n ), at intervals of 0.1 each, produces a 900 × 490 solution matrix containing all possible P aw values capable of satisfying Eq. 2 for given a set of (t k ), F aw (t k ), and PEEP a measurements made at time t k during insufflation. Figure 1 shows a schematic of the proposed numerical method of solution. In this example, a solution matrix was generated for ΔV(t k ) = 300 mL, F aw (t k ) = 32 L·min −1 , and PEEP a = 5 cmH 2 O and plotted as a three-dimensional surface in a Cartesian (C rs , R rs , P aw ) system. According to the above reasoning, the solution of Eq. 2, in terms of C rs and R rs , resides on a point on that surface. Further insight is gained by noting that the solution must lie along a surface path traced by the measured P aw (t k ) at time t k . This is shown in Fig. 1 as path A, where P aw (t k ) = 27 cmH 2 O and point 'a' symbolizes the yet unknown solution of Eq. 2. Fig. 1 Schematic of the numerical method used to solve the respiratory system equation of motion for static compliance (C rs ) and airway resistance (R rs ). In this example, the solution matrix was developed for ΔV(t k ) = 300 mL, F aw (t k ) = 32 L·min −1 , and PEEP a = 5 cmH 2 O and shown graphically as a three-dimensional surface bounded by C rs values ranging from 10 to 50 mL·cmH 2 O −1 and R rs from 0 to 20 cmH 2 O·s·L −1 . This surface encompasses all possible combinations of P aw , C rs and R rs capable of satisfying Eq. 2 for a given set of ΔV, F aw , and PEEP a measurements made at time t k during insufflation. P aw , also measured at t k and equal in this example to 27  Projecting path (A) onto the C rs -R rs plane generates a two-dimensional function, shown as path B, that restricts all possible combinations of C rs and R rs able to satisfy Eq. 2 for the set of measurements taken at time t k . The C rs -R rs function is developed numerically by noting the values R x , C y associated with those matrix elements where P aw (R x , C y ) = measured P aw (t k ).
It only remains to identify the location of solution point 'b' on the C rs -R rs plane. This is accomplished by generating a family of C rs -R rs functions, one for each set of P aw (t k ), ΔV(t k ), F aw (t k ), and PEEP a values measured at sequential times t k during insufflation. Since the one-compartment model of Eq. 2 assumes constant C rs and R rs , it follows that all generated C rs -R rs functions must pass through, and therefore intersect, at a point that defines C rs and R rs for the breath in question.
It is known that C rs and R rs vary early in inspiration as unstable alveoli open and conducting airways distend [10]. However, as lung volume increases past a lower inflection point (LIP) C rs achieves steady state until reaching an upper inflection point (UIP) where over-distention might occur. Defining the LIP and UIP by their respective lung volumes as ΔV LIP and ΔV UIP , it is reasonable to expect all C rs -R rs functions generated for insufflation lung volumes ΔV LIP < ΔV(t) < ΔV UIP to intersect at the solution point 'b' uniquely defining C rs and R rs for that breath. Figure 2 shows a family of C rs -R rs functions (n = 14) generated during a single breath's insufflation at sequential 32-ms intervals, past ΔV LIP > 200 mL. The intersection of these functions defines the values for C rs = 32.8 mL·cmH 2 O −1 and R rs = 23.8 cmH 2 O·s·L −1 . The inset graph illustrates the slight uncertainty associated in determining the intersection of the C rs -R rs functions, likely the result of random variations in measurement or small changes in C rs and R rs occurring during the insufflation. Accordingly, the point of intersection is best defined by the smallest standard deviation (σ) of all C rs values measured at each R rs increment along the R rs axis.

Fig. 2
A family of C rs and R rs functions (n = 14). Each function was generated at different times (t k ) measured sequentially at 32 ms during a single insufflation. The intersection of these functions defines C rs = 32.8 mL·cmH 2 O −1 and R rs = 23.8 cmH 2 O·s· L −1 for the breath. Shown in the inset graph is the uncertainty associated with the intersection point, likely the result of measurement limitations or minute alterations in C rs and R rs during insufflation. Accordingly, the point of intersection is best defined by the smallest standard deviation (σ) of all C rs values measured at each R rs increment along the R rs axis Gutierrez Intensive Care Medicine Experimental (2022) 10:55

Methods
The accuracy of the dynamic method was tested by comparing paired C rs and R rs values predicted by the dynamic method and the static method (used here as the 'gold standard') for the same breath.

Experimental validation
Validation was performed experimentally with a Maquet 190 one-liter test lung (Getinge, Solna, Sweden) using VC ventilation with a 0.5-s inspiratory hold. The test lung was attached to a Servo s ventilator (Getinge, Solna, Sweden) and ventilated at a respiratory rate of 15 bpm with V tidal of 300 mL or 400 mL. PEEP levels of 0, 5 and 10 cmH 2 O were applied sequentially at each V tidal . An in-house built data acquisition monitor was used to sample F aw and P aw signals from the ventilator data-port at 32.25 Hz and compile successive epochs of 4096 points, each lasting 131.1 s. Five epochs were obtained at each V tidal -PEEP combination. Data were analyzed in situ with the monitor's Raspberry Pi 3B processor programmed (Python 3.8) to calculate C rs and R rs for each breath by the dynamic method. C rs and R rs were also determined manually by the static method for 10 breaths in each epoch using data from the P aw and F aw signals. Average epoch values for C rs and R rs computed with either method were compared at each V tidal -PEEP combination.

Clinical validation
The dynamic method also was validated with clinical data using F aw and P aw signals obtained in a prior study of mechanically ventilated patients performed in 2011-2012 at The George Washington University Hospital Intensive Care Unit (IRB No. 110910) [11]. The database (Additional file 1: Section S1) contains information from 176 patients with acute respiratory failure enrolled within 24 h of intubation and monitored during their entire time on ventilatory support. It contains deidentified demographic information and F aw and P aw signals sampled at 32.25 Hz by the ventilator (Servo I or Servo S ventilators, Getinge, Solna, Sweden). The signals were saved as contiguous time-windows or epochs, each lasting 131.1 s and containing 4096 samples of each signal.

Epoch selection
Software was written (Python 3.8) to search the database for epochs on VC ventilation. The respiratory rate variability (RRV) for each identified epoch was used to determine the degree of active respiratory muscle activity. RRV was determined from the frequency spectrum of the expiratory flow signal as previously described [12] using the fast Fourier transform (FFT) algorithm [13]. RRV was defined as 100 -H1/ DC %, where H1 is the amplitude of the spectrum's first harmonic and DC that of the zero-frequency component. Epochs with RRV < 55% were assumed to have negligible respiratory muscle activity (P mus = 0) and were chosen for the study. This RRV value corresponds to those noted in normal individuals during quiet breathing in stages N2 and N3 of sleep [14].

Breath selection
Within each selected epoch, the software further identified breaths displaying an endinspiratory hold and absent voluntary respiratory effort. These breaths allowed for the reliable measurements of static compliance and airway resistance using standard calculations for comparison with those predicted by the dynamic method. The following criteria was used to choose breaths for analysis: (1) a discernible end-inspiratory hold > 0.25 s with mean plateau airway flow < 1 L·min −1 ; (2) ventilator-triggered (PEEP a -minimal P aw < 0.3 cmH 2 O); (3) full volume breaths (ΔV(t) ≥ 300 mL with insufflation time (Ti) > 0.8 s); (4) absent PEEP i (end-exhalation (EE) F aw < 3 L·min −1 [15] and P aw (t 0 )-EE P aw < 2 cmH 2 O) [16]; and (5) no air circuit leak (inspiredexpired V tidal <|30 mL|). Excluded were breaths with ΔV(t) ≥ 740 mL to avoid exceeding the UIP [17] (see Additional file 1: Section S2 for breath exclusion example). Once a breath was deemed adequate for analysis, the software calculated C rs and R rs by both the dynamic and static methods.
Results from the dynamic and static methods were compared with Pearson's linear regression R 2 and Bland-Altman analysis [18] for bias ± 95% limits of agreement (LOA). Since some patients had substantially more epochs meeting study criteria than others, the methods were compared in aggregate by combining data from all epochs and individually for each study patient. Unless otherwise specified, data are shown as median and interquartile range. The Mann-Whitney test was used to determine significant differences between independent samples. All reported p values are two-sided with p < 0.05 considered significant.

Experimental validation
Analysis of 880 breaths from 30 epochs resulted in nearly identical values for C rs and R rs calculated by the static and dynamic methods for all tested combinations of V tidal and PEEP (  Table S1, Additional file 2, Additional file 3, Additional file 4).
Of the 33,371 identified epochs, 3174 (9.5%) contained breaths displaying end-inspiratory holds. The ventilatory parameters associated with these epochs were compatible with those of quiet, passive ventilation with a low RR = 11 [11,14] bpm and RRV = 45 [40,46] % (Additional file 1: Table S2). The 3174 chosen epochs encompassed 87,021 individual breaths with 78,371 (90.1%) considered adequate for analysis of static compliance and airway resistance using standard calculations for comparison with those predicted by the dynamic method. The median number of breaths in these epochs was 24 [21,30].
Individual patient analysis.

Discussion
Increases in computing power [19] allow for the application of powerful analytical techniques to monitor patients on ventilatory support. The present study describes an algorithm capable of providing breath-by-breath measures of C rs and R rs without the need for an end-inspiratory pause in patients on VC ventilation. This technique may in turn allow for the continuous monitoring of other parameters, such as the driving pressure, a proven indicator of ventilator associated lung injury [20].
The dynamic method used to determine static C rs and R rs is based on a novel numerical solution of the equation of motion of the respiratory system. This equation depicts the behavior of respiratory mechanics in normal individuals and has been applied successfully to ventilated patients with respiratory failure [21]. In the form used here, the equation of motion ignores the inspired gas inertia and the resistance to energy transfer by visco-elastic lung tissue, whereas both terms may be quantitatively significant under extreme ventilatory conditions, they are likely inconsequential under the studied conditions [22]. It should be noted that the one-compartment model of Otis et al. [7] does not allow for the partitioning of respiratory system mechanics. On the other hand, a similar numerical approach may be considered when solving a more complex model of the respiratory system, one that accounts for both lung and chest wall compliances.
Method validation was done with matching pairs of C rs and R rs calculated by the static and the dynamic methods. Experimental method validation yielded nearly identical C rs and R rs values when tested with a test lung ventilated using different V tidal and applied PEEP levels. Since C rs and R rs were more or less fixed for the test lung, the dynamic method also was validated using airway signals from a previous study on mechanically ventilated patients. The use of clinical data yielded a more realistic assessment of the dynamic method, allowing for method comparison at C rs values ranging from 20 to 70 mL·cmH 2 O −1 and from 10 to 32 cmH 2 O·s· L −1 for R rs . These are ranges similar to those encountered in clinical practice. Software was written to identify breaths meeting strict morphologic criteria that included a discernible plateau pressure and negligible P mus or PEEP i . This resulted in the evaluation of a massive number of individual breaths (78,371) contained in the 3174 identified epochs. The software calculated paired C rs and R rs values by the static and dynamic methods in all identified breaths enclosed within each 131.1-s-long epoch, reporting the epoch's average for comparison. The use of epochs was dictated both by the format initially used to store the data and by the ability to assess respiratory muscle activity indirectly by spectral analysis. The cohort was composed mainly of highly sedated patients transferred from the Emergency Department and ventilated with end-inspiratory holds that were not immediately detected by the ICU team. Although the data were collected several years ago, neither the passage of time nor changes in ICU care should have influenced the results presented nor adversely altered the fidelity of the stored airway signals.
To provide for a balanced assessment of the data, analysis was performed in aggregate form and also individually for each patient in the study. Whereas aggregate analysis biased the results in favor of patients with many analyzed epochs, individual analysis amplified the effect of patients with fewer epochs. Regardless of comparison strategy, however, both methods produced nearly identical C rs and R rs values with negligible bias and exceedingly small LOA.
Although method bias was minimal for both C rs and R rs , the possibility should be acknowledged of introducing a systematic error by the software when calculating the "gold standards" C rs and R rs by the static method. The cessation of gas flow during the end-inspiratory hold produces a rapid decline in P aw from P peak to P 1 , followed by a slow decay to a plateau P 2 [23]. The timing of the end-inspiratory hold (t hold ) could be an important source of measurement error since a short t hold may affect P 1 by the persistence of airflow during inspiratory valve closure or P 2 by prematurely shortening the decay of P aw . Conversely, a long t hold may allow voluntary respiratory muscle activity to occur, also distorting P 2 . All breaths in the study were ventilator triggered with no evidence of spontaneous respiratory muscle activity throughout the length of the breath, including the end-inspiratory hold portion. For the cohort, t hold was 0.   [24] and placing P 2 firmly on the flat portion of the plateau, as evidenced by the small decline in P aw (< 1.0 cmH 2 O) predicted by decreasing exponentials fitted to the data (R 2 = 0.96) and extrapolated from 0.4 to 1.0 s (Additional file 1: Section S4). Several assumptions were made in the development of the dynamic method, among them the constancy of C rs and R rs during insufflation. This basic tenet of the onedimensional model of Otis et al. [7] is unlikely to hold true during the early stages of inspiration where the volume signal is curvilinear [23]. Past a certain inflation volume, defined here as ΔV LIP , C rs becomes constant and remains so over the rest of the tidal range [25]. The dynamic method was therefore applied to insufflation lung volumes > 200 mL, a ΔV LIP chosen to match those reported in ARDS patients [26,27]. This is probably a conservative estimate since no patient in the study met the Berlin definition for ARDS [28] with half the cohort having normal chest radiographs. Moreover, all patients were ventilated with PEEP a = 5 cmH 2 O, likely resulting in initial lung volumes in the region of constant C rs . It is possible, however, that small variations in C rs and R rs during the studied insufflation volumes resulted in the slight uncertainty noted in determining the intersection of the C rs -R rs functions.
The assumption of absent patient inspiratory effort during insufflation (P mus = 0) cannot be independently verified since esophageal balloon catheters were not used in the original study. The validity of this assumption rests on: (1) the use of RRV < 55% as an inclusion criterion, a value noted in heavily sedated ventilated patients [29] and normal individuals during stages of deep sleep; (2) all analyzed breaths were ventilator-initiated; and 3) a cohort of 50 epochs selected randomly from the sample population was characterized by a regular breathing pattern, low respiratory rate (11 [11,14] bpm) and no signal distortion (see Additional file 1: Section S7 and Table 2e).
The assumption of absent PEEP i also could not be independently verified, but care was taken to include in the analysis only breaths displaying minimal differences between its onset and the prior breath's end-exhalation F aw and P aw . In addition, (1) no patient in the study was diagnosed with obstructive lung disease; (2) the exhalation time for the cohort allowed ample time for expiration (3.2 ± 0.7 s); and (3) tachypnea (RR > 20 bpm) was absent in all chosen epochs.
The dynamic method is unlikely to perform well under conditions of persistent asynchronous breathing or in the presence of significant respiratory muscle effort. It is also not amenable for bedside use or with ventilators lacking airway signal sampling. Conversely, when used in conjunction with a computer connected to the ventilator's data-port, the dynamic method may provide accurate ongoing measurements of C rs and R rs under most clinical conditions encountered during the provision of volume-controlled mechanical ventilation.
Although the present study was not intended as a methodological comparison, the dynamic method appears to perform as well or better than either the LSF or the τ e methods (Additional file 1: Table S3). Unlike these empirical models, the dynamic method represents a deterministic approach to the solution of the equation of motion. As such, it may be applicable to ventilatory modes other than VC and provide insight into the relationship of respiratory system mechanics to other ventilatory variables, such as plateau pressure, respiratory muscle effort and intrinsic PEEP. These, and other issues related to the application of the dynamic method await further confirmation by prospective studies. (2022) 10:55 Take-home message A novel numerical method to calculate static compliance and airway resistance of the respiratory system during ventilatory support is developed and validated.